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Abstract. The cosmic-ray streaming instability creates strong magnetohydrodynamic turbulence in the precursor 
^ ' of a SN shock. The level of turbulence determines the maximum energy of cosmic-ray particles accelerated by 

1/^ , the diffusive shock acceleration mechanism. The present consideration continues our work Ptuskin & Zirakashvili 

■ (2003). It is assumed that the Kolmogorov type nonlinear wave interactions together with the ion-neutral collisions 

^— ' restrict the amplitude of random magnetic field. As a result, the maximum energy of accelerated particles strongly 

depends on the age of a SNR. The average spectrum of cosmic rays injected in the interstellar medium over the 
course of adiabatic SNR evolution (the Sedov stage) is approximately Q{p)p^ (x at energies larger than 10 — 30 
Gev/nucleon and with the maximum particle energy that is close to the knee position in cosmic ray spectrum 
observed at ~ 4 x 10^^ eV. At earlier stage of SNR evolution - the ejecta dominated stage described by the 
Chevalier-Nadyozhin solution, the particles are accelerated to higher energies and have rather steep power-law 
distribution on momentum. These results suggest that the knee may mark the transition from the ejecta-dominated 
^ ■ to the adiabatic evolution of SNR shocks which accelerate cosmic rays. 

, Key words. ISM: supernova remnants - ISM: cosmic rays -acceleration of particles - shock waves -turbulence 

gamma rays: theory 
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^ 1. Introduction the gas with density uq cm~^ and the interstellar magnetic 

field Bq — 5 fiG is considered. The found value of E'max 

... The diffusive shock acceleration is considered as the main ^j^^^ somewhat less than the energy of the "knee" , 

mechanism of acceleration of galactic cosmic rays. The ^^^^j^ ^^^^1 ^^^^-^ spectrum observed at 

dependent on energy diffusion coefficient D{E) determines ^ 4 ^ 10^^ eV 

the maximum energy that particles can gain in the process 

of acceleration by the shock moving through the turbulent Analyzing the early stage of SNR evolution when the 

interstellar medium. The condition of efficient acceleration shock velocity is high, u^h ~ W km s"\ Bell & Lucek 

is D{E) < >fMshi?sh, where R^b is the radius and u^h is dZOOlj found that the cosmic-ray streaming instabihty 

the velocity of spherical shock, the constant >f - 0.1, see in the shock precursor can be so strong that the ampli- 

Drury et al. iPMjl . Malkov & Drury gOOlJ for review, fied field SB > 100 ^iG far exceeds the interstellar value 

The Bohm value of the diffusion coefficient Dg = vr^/i Bq. The maximum particle energy increases accordingly. 

{v is the particle velocity, and is the particle Larmor The cosmic-ray streaming instability is less efScient as 

radius) that is a lower bound of the diffusion along the the shock velocity decreases with time and the nonlin- 

average magnetic field, gives the maximum particle energy wave interactions reduce th e level of turbule nce a t 

-E,„ax 2 X 10"Z(£'5i/no) eV at the time of transition the late Sedov stage (Volk et al. HHEHI Fedorenko [111111) . 

from the ejecta dominated stage to the stage of adiabatic This leads to the fast diffusion and to the corresponding 

evolution of SNR (the particle charge is Ze). Here the SN decrease of E^^^. The effect is aggravated by the possi- 

burst with the kinetic energy of ejecta £ = £51 lO^^ erg in ble wave dan iping on the ion-neutral colhsions (Bell 119781 

Drury et al. 'TSUB"). The acceleration of cosmic rays and 

Send offprint requests to: V.S. Ptuskin their streaming instability in a wide range of shock veloc- 
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ities was considered in our paper Ptuskin & Zirakashvili 
H2003|l (Paper I). The analytical expressions for cosmic 
ray diffusion coefficient and for the instability growth rate 
were generalized to the case of arbitrary strong random 
magnetic field, 6B ^ Bq. The rate of nonlinear wave in- 
teractions was assumed to correspond to the Kolmogorov 
nonlinearity of magnetohydrodynamic waves. The coUi- 
sional dissipation was also taken into account. The maxi- 
mum energy of accelerated particles was determined as a 
function of shock velocity and thus as a function of SNR 
age. The maximum energy of the particle with charge Ze 
can be as high as lO^^Z eV in some very young SNRs 
and falls down to about 10^°Z eV at the end of adiabatic 
(Sedov) stage. The widely accepted estimate of cosmic ray 
diffusion coefficient at the strong shock that corresponds 
to the Bohm diffusion value calculated for the interstellar 
magnetic field strength turns out to be not correct. This 
result may explain why the SNRs with the age more than 
a few thousand years are not prominent sources of very 
high energy gamma-rays (Buckley et al llH98l Aharonian 
et al I2()()2|l . The presence of strongly amplified random 
magnetic field in young SNRs is evidently supported by 
the interpretation of data on synchrotron X-ray emission 
from young SNRs, see e.g. Vink H2003|l for review. 

The main objective of the present work is the calcu- 
lation of the average spectrum of cosmic rays ejected in 
the interstellar medium by a SNR in a course of its evolu- 
tion. Some necessary results of Paper I are presented in the 
next Section 2, the evolution of SNR shocks is discussed in 
Section 3, the average cosmic-ray source spectrum is calcu- 
lated in Section 4 followed by the discussion in Section 5, 
the conclusion is given in Section 6. Appendix A describes 
the thin shell approximation used in our calculations. 

2. Maximum Energy of Accelerated Particles 

In the test particle approximation, the distribution of ac- 
celerated particles in momentum for high Mach number 
shocks has the canonical form /(p) ~ fKrvmskv ll977l 
Bell 1978J . In the case of efficient acceleration, the action 
of cosmic ray pressure on the shock structure causes non- 
linear modification of the shock that changes the shape 
of particle spectrum making it flatter at ultra relativis- 
tic energies (Eichler 1984, Berezhko et al. 119961 Malkov 
& Drurv I^Olfl . Because of this effect, we assume that 
the distribution of ultrarelativistic particles at the shock 
is of the form fo{p) ^ p~^^°- where < a < 0.5, and 
the value a — 0.3 is used in the numerical estimates 
below. The normalization of function f(p) is such that 
the integral N = An J dpp^f{p) gives the number den- 
sity of cosmic rays. The differential cosmic ray intensity is 
I{E) — f{p)p'^. We assume that the cosmic ray pressure at 
the shock is some fraction ^cr < 1 of the upstream momen- 
tum flux entering the shock front, so that Per = ^cvPUsh 
and the equation for the distribution function of relativis- 
tic accelerated particles at the shock is 

f( ,^ ^ 3gcrPMsh-H^bmax(0 -P) 

4^c(mc)-^(p,„ax)p4-» ' ^ ^ 



where Pmax is the maximum momentum of accelerated 
particles, H{p) is the step function, and (^(pmax) — 



Jo 



^/mc dyy" 



The approximation of the last integral 



yiT^ 
at Pmax S> mc is (p{p) ~ a^^lp/mc)"' — a'~^{l + a)^^. The 
value of ^cr ~ 0.5 and the total compression ratio at the 
shock close to 7 were found in the numerical simulations 
of strongly modified SN shocks by Berezhko et al. H1996fl . 
Here and below we mainly consider protons as the most 
abundant cosmic ray component. For ions with charge Z , 
the equations should be written as functions of p/Z in- 
stead of p. In particular, the nuclei with charge Z reach 
the maximum momentum a factor of Z larger than pro- 
tons. We use the notation m for the proton mass. The 
acceleration in old SNRs (^ ^ 3 x 10* — 10^ yr) when 
Pmax/mc < 10 are not considered in the present paper 
because Eq. (1) is not applied at low Mach numbers, see 
Paper I for detail. [Using the test particle approximation 
for not modified shock, Drury et al. H2003fl found that the 
spectrum of accelerated particles is somewhat steeper if 
the diffusion coefficient is increasing with time compared 
to the case of constant D. This effect is not included in 
our consideration.] 

The following steady-state equation determines the en- 
ergy density W{k) {k is the wave number) of the mag- 
netohydrodynamic turbulence amplified by the streaming 
instability in the cosmic-ray precursor upstream of the su- 
pernova shock: 



.VW^(fc) = 2(r„ - Pi - r„i)iy(A:). 



(2) 



Here the l.h.s. describes the advection of turbulence by 
highly supersonic gas flow. The terms on the r.h.s. of the 
equation describe respectively the wave amplification by 
cosmic ray streaming, the linear damping of waves in back- 
ground plasma, and the nonlinear wave-wave interactions 
that may limit the amplitude of turbulence. The equation 
for wave growth rate at the shock 



r„(A:) = 



cr f^sh 



(1 



(l-a)/2 



cVa(/3(Pmax)?'°o 



(3) 



was suggested in Paper I as the generalization of equation 
derived for the case of weak random field (Berezinskii et al. 
11990^. Here T4 = Bo/^/Anp is the Alfven velocity (p is the 
gas density), A = SB/Bq is the dimensionless wave ampli- 
tude, and Tgo = mc^/eBa. The ion-neutral and electron- 
ion collisions usually determine the linear damping pro- 
cesses in the thermal space plasma. The Kolmogorov-type 
nonlinearity with a simplified expression 



{2CKy^^^V^kA{> k) w 0.05Kfc^(> k) 



(4) 



at Ck — 3.6 (as it follows from the numerical simula- 
tions by Verma et al. I1996|l was used in Paper I. The 
wave-particle interaction is of resonant character and the 
resonance condition is k^c^fg — -\/l+^4^, where the 
Larmor radius is defined through the regular field Bq 
Tg = pc/ZeBo , and Atot is the total amplitude of ran- 
dom field. The particle scattering leads to the cosmic- 
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ray difFusion with the diffusion coefficient D = (1 + 
Al^.Y^-'vT^ [3A2(>fc,,,)]-\ 

Eq. (2) ailows finding the foilowing approximate equa- 
tion for the dimensionless amphtude of the total random 
magnetic field produced by the cosmic-ray streaming in- 
stability at the shock (Paper I): 



2v{\ + Al,,) (2Ck)3/2 

Ccr(Q)CcrMsh 

cVa<y5(Pmax) (Pmax/»71C)-" yr 



To study the effect of nonlinear interactions, the term with 
linear damping in Eq.(2) was omitted in Eq.(5); Ccr(a) = 
27[4(5 — a) (2 — a)]~^. The maximum particle momentum 
satisfies the equation 



(6) 



In the high velocity limit, when itgh S> 
4aC„^erc[9(2C'K)3/2]-i and » W,[2aC„S,,,r\ 

the advection term dominates over the nonlinear dissipa- 
tion term in the l.h.s. of Eq. (5) and the wave amplitude is 
large, ^tot 3> 1. The maximum momentum of accelerated 
particles and the amplified magnetic field are given then 
by the approximate equations 



Pn 



mc 
and 



A, 



2Msh 

3K 



(7) 



(8) 



The cosmic ray diffusion coefficient depends on particle 
mor radius as I? cx vrl^'^ at p < Pmax in this case. 
In the low velocity limit, when Ugh ^ 



W^c' (7raC„(2CK)^ecr)- 

K?c(aC„(2CK)3/2^er)"' 



1/5 
1/3 



and 



the nonlinear dissipa- 
tion term dominates over the advection term in the l.h.s. 
of Eq. (5) and the wave amplitude is small, Atot 1- 
The maximum momentum of accelerated particles and 
the amplified magnetic field are given then by the 
approximate equations 



Pn 



mc 
and 



(9) 



(10) 



assuming that this value of amplified field exceeds the 
value of random interstellar magnetic field. The cosmic 
ray diffusion coefficient depends on particle Larmor ra- 
dius as Z? cx wrg~^° at p < Pmax- (Notice the misprint in 
the numerical coefficient in the first equality of Eq. (19) 
in Paper I that is analogous to the present Eq. (9).) 



Fig.l illustrates the results of calculations of Pmax at 
the Sedov stage of SNR evolution at £ = 10^^ erg in the 
warm interstellar gas with the temperature T — 8 x 10^ 
K, the average density no = 0.4 cm""^ that includes small 
interstellar clouds, the intercloud density n = 0.1 cm~^, 
the number density of ions = 0.03 cm""^, the interstellar 
magnetic field value Bq = 5 /iG, see Paper I. The time 
dependence of the shock radius and the shock velocity are 
given by the following equations (the Sedov solution, see 
e.g. Ostriker & McKee HM^ : 



(5) i?sh = 4.3(£5i/no)'/'C 



pc. 



= 1.7 X 10^ (£5l/"0)'/'tKyr' 



(11) 



where we assume that the ultrarelativistic gas of cosmic 
rays mainly determines the pressure behind the shock. The 
value K — 0.0 A was assumed in the calculations in Fig. 1. 
Three solid lines correspond to the three cases of wave dis- 
sipation considered separately: the nonlinear wave inter- 
actions; the damping by ion-neutral collisions at constant 
gas density; the damping by ion-neutral collisions when 
the diffuse neutral gas restores its density after complete 
ionization by the radiation from the SN burst. For the last 
two curves, the dissipation of wave due to the ion-neutral 
collisions with damping rate 



1 + A 



(1 + -)^ 



(12) 



was taken into account whereas the term T^i that de- 
scribes the nonlinear dissipation was omitted. Here fin = 
riH {vth<y) « 8.4 X 10-9(T/104 K)° *{nii/1 cm^^) g-i for 
the temperature T ~ 10^ — 10^ K is the frequency of ion- 
neutral collisions with the cross section a averaged over 
velocity distribution of thermal particles, njj is the num- 
ber density of neutral hydrogen. (Notice some correction 
of the "coUisional damping" curve in Fig.l comparing to 
corresponding Figure 2 in paper I.) The maximum energy 
of protons accelerated by SN shocks at the early Sedov 
stage is close to 3 x 10^** eV that exceeds the Bohm limit 
calculated for the interstellar magnetic field value by one 
order of magnitude. The maximum energy decreases to 
about 10^" eV at the end of the Sedov stage that is much 
less than the Bohm limit calculated for the interstellar 
magnetic field value. In particular, the particle energy is 
less than lO^'^ eV at i > 3 x 10'^ yr and this may explain the 
absence of a TeV 7-ray signal from many SNRs (Buckley 
et al.il998| Aharonian et al. 2002) where the gamma-rays 
could in principle be produced through tt" decays if suffi- 
ciently energetic cosmic rays were present. 

With the extreme choice of parameters of the fly- 
ing apart young SNR envelope, it was found (Bell & 
Lucek EE 



m Paper I) that the maximum particle energy 
may reach the ultra high energies. The estimate of the 
highest particle energy according to Paper I is ii^max ~ 
2 X 10i^Z(wsh/3 X W km s-^)'^{K/OA)^„M^/^n^^^ eV at 
the end of a free expansion stage which precedes the Sedov 
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no wave damping 




SNR age 
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recombination 
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Fig. 1. The maximum momentum of accelerated protons pmax in units mc as a function of shock velocity Ush at the 
Sedov stage of supernova remnant evolution in warm interstellar gas. Three solid lines correspond to three cases of wave 
dissipation considered separately: nonlinear wave interactions; damping by ion-neutral collisions at constant density 
of neutral atoms; damping by ion-neutral collisions when the diffuse neutral gas restores its density after complete 
ionization by the radiation from supernova explosion. The dashed line presents the age of a supernova remnant t 
(plotted on the right ordinate) as a function of shock velocity. The dotted line shows the Bohm limit on maximum 
particle momentum calculated for the interstellar magnetic field strength. The dash-dot line gives the maximum 
particle momentum when the wave dissipation is not taken into account. 



stage (here Mej is the mass of ejecta in solar masses). We 
shall see below that this promising estimate is in some 
sense devaluated by the results of calculations of particle 
flux - the flux turns out to be low at the highest energies 
which can be achieved in the process of acceleration. 



3. Evolution of SNR Shocks 

The typical source of galactic cosmic rays is most probably 
associated with the core collapse supernova, Type II SNe, 
that is the final stage of evolution for stars more massive 
than about 8 solar masses while on the main sequence. 
The massive star before the explosion goes through the 
sequence 0-star stage, Red Super Giant star stage, and 



through the Wolf-Rayet stage for the most massive pro- 
genitors (> 2OM0) that give the rare in occurrence Type 
Ib/c SNe. The fast wind of a massive progenitor star on 
the main sequence produces a big bubble of hot rarefied 
gas with the temperature about 10^ K in the surrounding 
interstellar medium, see Weaver et al. (|1977f) . Lozinskaya 
(UnnUl- The typical Type II SN goes through the Red 
Super Giant phase before the explosion and this process 
is accompanied by the flow of a low-velocity dense wind. 
Thus, immediately after the supernova burst, the shock 
propagates through the wind of a Red Super Giant star 
then through the hot bubble and finally it enters the in- 
terstellar medium. Our calculations will be done for the 
ejecta mass M^. = IMq (the solar mass). The spherically 
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symmetric distribution of gas density in the stellar wind 
is Tiw = M/(47rTOaWw?'^), where M — 10~^Af_5 (solar 
mass)/yr is the mass loss rate, ma = 1.4m is the mean 
interstellar atom mass per hydrogen nucleus, the wind 
velocity Mw = ^O^u^fi cm/s. The stellar wind magnetic 
field has the shape of the Parker spiral similar to the case 
of interplanetary magnetic field (Parker, 1958). At rela- 
tively large distances from the surface of the star that are 
of interest here, the magnetic field has predominately az- 
imuthal structure and its value is Bq = B^,r'^^} am 9 /{uwr) 
where is the surface magnetic field strength at the star 
radius r*, J7 is the angular velocity of star rotation, and 
9 is the polar angle. Hence Bo{r)r — 2 x 10^'^M~gSin^ 
Gxcm at = 1 G, = 3 X lO^^ cm, 17 = 3 x 10"* s^^ 
that gives i?o ~ 6 /iG at the distance r = 1 pc from the 
star. 

Below we shall also use the following set of parameters 
of the medium surrounding the Type II SN: the radius 
of spherical Red Super Giant wind R„ — 2 pc, the star 
mass loss M_5 = 1, and the wind velocity m^.q = 1- The 
radius of the spherical bubble of hot gas Rb = 60 pc , 
the gas density in the bubble rib = 1-5 x 10~^ cm~^, the 
magnetic field there i?b = 5 /xG. The gas density in the 
undisturbed interstellar medium around the bubble is as- 
sumed to be equal to ng = 1 cm~'^ (physically, the value of 
no determines rib, see Weaver et al. 1977). The hot bubble 
is separated by the dense thin shell from the interstellar 
gas. The accepted parameters are close to those assumed 
by Berezhko & Volk (20QQ| in their analysis of gamma- 
ray production in SNRs. The lengthy discussion and the 
additional references can be found there. 

The considerable fraction of cosmic rays is probably 
accelerated in Type la SNe (their explosion rate in the 
Galaxy is about 1/4 of supernovae Type II rate). These 
supernovae are caused by the thermonuclear explosions of 
compact white dwarfs following mass accretion. The char- 
acteristic mass of a progenitor star and the mass of ejecta 
are 1.4 solar mass. The progenitor stars do not appear 
to have observable amount of mass loss nor do they emit 
ionizing radiation that could essentially modify the am- 
bient medium around the star. We assume that the SNR 
shock goes through the uniform weakly ionized interstel- 
lar medium with density 1 cm"^, and the magnetic field 
7 AiG. 

The two asymptotic regimes of the propagation of SNR 
shock - the ejecta dominated stage and the adiabatic stage 
- are instructive to consider. 

The adiabatic regime was mentioned earlier, see the 
Sedov solution (11) for the shock moving in the gas with 
constant density, and it refers to the stage of SNR evo- 
lution when the mass of swept-up gas significantly ex- 
ceeds the mass of ejecta. This condition is fulfilled in the 
case of the medium with constant density at Rsh > Rq = 
(3Mej/47rmano)^/^ = l.9{M^j/MQno)^/^ pc, to > R/uq « 
190nQ ^^'^ yr, where uq ~ 10^ cm/s is the ejecta initial 
velocity. The adiabatic regime for the SNR shock mov- 
ing through the progenitor star wind is described by the 



equations (at Ush ^ "w): 

/ c \ 1/3 

D ~ n '^5lUwfi \ ,2/3 

= 5.2 X 10^ (y^Yrf) ^ ^"^''^ ^^^^ 

see Ostriker & McKee (1998). As in Eq. (11), we assume 
that the ultrarelativistic gas of cosmic rays mainly deter- 
mines the pressure behind the shock. Eq. (13) is vahd 
when the mass of swept-up gas is relatively large and 

i?sh > i?0 = MejUw/M W l(Mej/MQ)Uw,6/M_5pC. 

The quantity pu^^R^^ = K£ is conserved for the con- 
sidered adiabatic shocks. The constant K « 0.16 for 
the solution (11), and K w 0.34 for the solution (13). 
In the general case of the power-law gas distribution 
p = po{r)r~'^ , s < 5, the adiabatic shock evolution is de- 
scribed by the equations Rsh = {vis)^ / Po)'^^t'^^ , and 
Ush — -^-^ivi-^)^ / Po)'^^t^^^ : where rj is constant at 
fixed s, that gives the general formula K = (l^-gyi ( the 
values of rj{s) were given by Ostriker & McKee. HMSjl . 

The ejecta dominated stage precedes the adiabatic one. 
As long as the mass of the ejecta is large compared to 
the swept-up mass, the blast wave is moving with rela- 
tively weak deceleration. At this stage shortly after the 
explosion, the structure of the flying apart envelope of the 
presupernova star is important for the shock evolution. 
Actually, the blast wave consists of two shocks, the for- 
ward shock and the reverse shock, with the contact dis- 
continuity surface between them. This surface separates 
the shocked wind or interstellar gas downstream of the 
forward shock from the shocked envelope gas that fills the 
downstream region of the reverse shock. The reverse shock 
lags behind the forward shock and enters the dense inter- 
nal part of the fiying apart star by the time of the begin- 
ning of Sedov stage. Though it can not be well justified ap- 
proximation for the very young SNRs, we ignore below the 
cosmic ray acceleration at the reverse shock compared to 
the forward shock and use the notation Rsh for the radius 
of forward shock (Berezinsky & Ptuskin (1989) consid- 
ered the cosmic ray acceleration by both shocks, see also 
Yoshida & Yanagita H20()l|l V The outer part of the star 
that freely expands after the SN explosion has a power 
law density profile ps cc r"*^, see e.g. Chevalier & Liang 
H1989|l . The value of k typically lies between 6 and 14. 
The value A: « 7 is characteristic of the SNe Type la, and 
A; « 10 is typical of the SNe Type II. The self similar solu- 
tion for the blast wave at the ejecta dominated stage was 
found by Chevaher 1(1^5^ and Nadyozhin (1981, 1985|. It 
was shown that at the age larger than about one week, the 
evolution of the shock at the ejecta dominated stage can 
be approximately described by the power-law dependence 
Rsh oc where the expansion parameter A = for the 
explosion in the uniform medium, and A = for the 
explosion in the wind of a presupernova star (for k > 5 
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ejecta). In particular, using the results of two mentioned 
above papers, one can obtain the following equations 



Rsh = 5.3 



(14) 



for the Type la SN explosion in the uniform interstellar 
medium at fc = 7; 

/ 7/2 5/2 \ "'"/^ 



5/2 



Ush — 6.6 X 10 



p7/2 »#5/2\ 1/8 

^' ^""-'^t^f t-^li'^m/s (15) 



for the Type II SN explosion in the wind of a presupernova 
star at fc = 10. 

Following the approach of Truelove & McKee (1999), 
one can describe the shock produced by the Type la 
SN using the continuous solution which coincides with 
the ejecta dominated equation (14) until the moment 
to = 260{M^i/lAMQ)^/^£~^^^^n~^^^ yr, and is given by 
the equations 



i?,h = 4.3(f5l/?^0)'/' tKyrX 



^ 0.06(Mej/MQ) 



5/6 



2/5 



cl/2„l/3, 
<--51 "o ''Kyr 



pc, 



(16) 



Ush = 1.7x 10^ (f5l/no)'^'iKyr'x 
V Szi Uq' <Kyr / 

at a later time t > to. It is evident from Eq. (16) that the 
adiabatic asymptotics (11) holds at t ^ to- 

The evolution of the Type II SN shock first follows the 
ejecta dominated solution (15) in a presupernova wind and 
then, while still moving in the wind, it enters the adiabatic 
regime at the distance r ~ 1 pc. The subsequent evolu- 
tion proceeds in the medium with a complicated structure 
described above for the Type II SN. The fairly accurate 
solution for the SNR evolution during this period can be 
obtain in the "thin-shell" approximation, e.g. Ostriker & 
McKee (1988J, Bisnovatyi-Kogan & Silich (1995). Using 
this approximation for the strong shock and assuming the 
spherically symmetric distribution of the circumstellar gas 
density p(r), we come to the following equations where the 
shock velocity Ush and the SNR age t are parameterized 
as functions of the shock radius Rsh (see Appendix for the 
derivation of these equations): 



Ush (-Rsh) 



7ad 



1 



12(7ad - l)£ 



(7ad + l)M2(i^3h)i?'^'^''^/^^'^ + '^ 



_Rsh 



drr V ^ad+i. 



M{r) 



1/2 



(17) 



t(-Rsh) 



dr 



where -fad is the adiabatic index (7ad = 4/3 if the pressure 
downstream of the shock is determined by the relativistic 
particles), M{R) — Mej + 47r drr'^ p{r) is the mass of 
the swept up gas. The self similar solution by Chevalier 
and Nadyozhin is not exphcitly reproduced by Eqs. (17). 
The solutions (15) and (17) are fitted together at the tran- 
sition from the ejecta dominated regime to the adiabatic 
regime (at r ~ 0.3 pc) in our numerical simulations of 
cosmic ray acceleration in the Type II SNRs described 
below. 

It is worth noting that the energy loss of SNR in a form 
of escaping cosmic rays is not taken into account in the 
solutions for shock evolution that were described in this 
Sections. In fact the shock evolution is only approximately 
adiabatic. 



4. Average Spectrum of Cosmic Rays Injected in 
the Interstellar medium 

At a given SNR age t, the cosmic rays are accelerated 
up to maximum momentum Pmax(0- Also, particles with 
P > Pmax(i) cannot be confined in the precursor of the 
shock even if they were accelerated earlier. Thus particles 
accelerated to the maximum energy escape from a SNR 
(see also Berezhko & Krymsky 1988). Let us estimate the 
flux of these run-away particles. We consider the simpli- 
fied approach for maximum energy of accelerated particles 
and take the dependence of diffusion on momentum in the 
following simplified form: 



D{p) = Do « RUsh, P ^ Pmax(i), 
D{p) = RUsh, P > Pmax{t)- 



(18) 



The spectrum of accelerated particles in this case has 
a very steep cut-off at p > Pmax (cf. Eq (1)) and the 
spectrum of run-away particles beyond Pmax can be ap- 
proximated by ^-function. To find the equation for these 
particles, let us integrate the equation for cosmic-ray dis- 
tribution function 



(19) 



on momentum p from Pmax to Pmax + Ap, where Ap << 
Pmax and larger then the width of run-away particle spec- 
trum. Denoting G = JP'^'^^'^^p fdp one obtains from Eq 
(19): 



dG 
'dt 



uVG 



dp„ 



at 



-/(Pmax-O) 



Vu 



Pmax/(Pmax-0).(20) 
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Since diffusion coefficient of run-away particles is large, 
the advection terms play no role in this equation and the 
last two terms can be considered as source of particles. The 
total source of run-away particles is given by the volume 
integral of this terms. As a result, the source spectrum of 
run-away particles has the form 



Q(P,t) = -<^(P-Pmax)x 



dpn 



dt 



^Pmax ) f{p 



0,r) 



(21) 



The integration here is performed over the domain where 
the integrand is negative. The integral 47r / dpp^q{p, t) has 
dimensions number of particles per unit time. 

Below we consider the case of spherically symmetric 
SN shock with linear velocity profile at r < Rgh'. 



u = I 1 - ^ ) Ush{t)r/Rsh{t), 



(22) 



where a is the total shock compression ratio. It includes 
a thermal subshock and a cosmic ray precursor. Linear 
profile of velocity (22) is a good approximation of Sedov's 
solution and it can be considered as a very approximate 
one at the ejecta-dominated stage. Since the shock is par- 
tially modified in the presence of cosmic rays, we should 
not assume any relation between the shock compression 
ratio a and the spectral index of accelerated particles 
4 — a (recall that 4 — a = 3c7/(it — 1) for not modified 
shocks). We accept the value cr = 7 in our calculations. 
The preshock at r > i?sh is created by the cosmic ray 
pressure gradient. Its width is small in comparison with 
the shock radius under the conditions given by Eq. (18) 
and the plane shock approximation can be used. Since the 
cosmic ray pressure dominates the gas pressure in the pre- 
cursor region, its gradient is proportional to the velocity 
gradient dPcr/dr = pushdu/dr, where p is the circum- 
stellar medium density. We also use an assumption that 
cosmic ray pressure at the shock is some fraction ^cr of 
the upstream momentum flux, see Eq. (1). Now assuming 
that /(pmax) is proportional to the cosmic ray pressure 
the expression (21) for the run-away particle source takes 
the form 



q{p, t) = Ati5{p - Pmax) 



1 



1 6 



R'^Ushpfoip)' 



f 

Jo 



r'^drfipn.a.^, 



r) 



dpn 



dt 



1 Ush 



^ (23) 



The first term in this exprc!ssion describes the particles 
which runs away from the shock front, and the second 
term describes the particles escaping from the shock inte- 
rior. In principle, the turbulence downstream the strong 
shock might be enhanced that would result in the small 
cosmic ray diffusion coefficient. In this case the particles do 
not run-away from the downstream and the second term 
in Eq. (23) should be omitted. If the turbulence down- 
stream is maintained by the same process of the cosmic 



ray streaming instability as in the upstream region, the 
downstream diffusion coefficient is comparable to the up- 
stream diffusion coefficient for particles with p ^ Pmax- 
We shall furhter assume that particles can run away both 
from upstream and downstream of the shock. The uncer- 
tainty of the efficiency of run-away process in the inner 
part of SNR does not qualitatively change the conclusion 
about the average source spectrum of cosmic rays calcu- 
lated later in this Section and shown in Fig. 2. 

The distribution function of particles with p ^ Pmax 
can be found using the solution of transport equation (19) 
at r < i?sh with the boundary condition f{p,r = Rsh, t) = 
fo by the method of characteristics. As the result 



q{p,t) = 47r5(p-pmax) 



ih-l-|:)i?W/o(p)+ 



dpn 



dt 



a 



4im Mm) 



(24) 



The expression in brackets in front of the integral in Eq. 
(24) should be positive that means that the particles lose 
energy adiabatically slower then the maximum energy de- 
creases. For the opposite sign, the adiabatic losses of par- 
ticles are faster then the decrease of maximum energy 
and the particles don't run away from downstream of the 
shock. They can run away at later time if at that time the 
decrease of maximum momentum will be faster. 

The average source power Q{p) of run- away cosmic 
rays per unit volume in the galactic disk is obtained by the 
integrating q with respect of t and by the averaging over 
many SN explosions: Q{p) = Ugn J"^-^^ dtq{p,t), where 
is the average frequency of SN explosions per unit volume 
of the galactic disk. Changing the variable of integration 
from t to Rsh {dRsh = Ushdt) one can derive the following 
equation: 



Qip) 



d ln(pnia3i) 



dln(fl) 



1 ^ dlnij) 



dln{R) 



dR' p{R')ul,^{R')R''- 



+ 



H [pn,UR')-p[§ 



P 



Pmax(-R') 



(a-l)(<7-l)/<7- 



(25) 



J R=Rm{p) 



Here Eq. (1) for fo with the approximation ip{p) w 
a~^{p/mc)°' — a~^{l + a)^^ is used, and the condition 
^cr = const is assumed. The function Rm{p) in Eq. (25) 
is defined by the equation Pmax(-Rsh = -Rm(p)) = P- If 
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the last equation has multiple solutions at given p, the 
summation on all these solutions should be performed in 
(25). The physical meaning of i?m(p) is that it is the value 
of shock radius when the maximum energy of accelerat- 
ing particles is equal to p. The second term in the r.h.s. 
of Eq. (25) should be omitted if the expression in round 
parenthesis in front of the integral is negative. 

Let us assume that the maximum momentum is a 
power law function of the shock radius, pmax oc R^^ , the 
particles are ultrarelativistic, p >> mc, and the compres- 
sion ratio is constant, a = const. The remarkable feature 
of Eq. (25) is then that the expression in square brack- 
ets does not depend on momentum at the adiabatic stage 
of SNR shock propagation in the medium with a power- 
law distribution of gas density, because pUgi-^Rf^^ = K£, 
K = const in this case, see Section 3. The average cosmic 
ray source power is now given by the simple equation 



particles in the high velocity limit (7) has the scaling 



Q{p) 



cp^ 



1 

3d 



1 



1 



+ 



1 - 



(t5 



(26) 



Here the factors 6, and (1 — should be positive. The 
first term in the square brackets describes the particles 
which run-away from the shock and the second term de- 
scribes the particles which run-away from the SNR in- 
terior. Consequently, while in the adiabatic regime, the 
SNR shock during its evolution produces the run-away 
particles with the universal power-law overall spectrum 
Q{p) oc p~^, whereas the instantaneous spectrum at the 
shock is more flat and not universal (sec Eq. (1)) and the 
instantaneous spectrum of run-away particles has a delta- 
function form (see Eq. (24)). 

The total source power of ultrarelativistic parti- 
cles calculated with the use of Eq. (26) is = 
4ttc Jdpp^Q{p) = CCcri'sn^ In(pmax2/Pmaxi), where C = 



12-KKa ( ^ 



l_ j_ _ ^ 



Pmax2 



and Pmaxi are the maximum momenta of accelerated par- 
ticles at the beginning and at the end of the adiabatic 
stage respectively, thus typically In(pmax2/Pmaxi) ~ 10. It 
leads to the estimate W « 0.5(^cr/0.5);/snf for the shock 
moving in the uniform interstellar medium. Hence the con- 
siderable part of the total available mechanical energy of 
SN explosion Vsn£ goes to cosmic rays at ^cr ^ 0.5. As it is 
well known, the source spectrum oc p~^ or somewhat more 
steep, and the efSciency of cosmic ray acceleration at the 
level 10 — 30% are needed to fit the cosmic ray data below 
the knee in the cosmic ray spectrum at about 4 x 10^^ eV 
in the empirical model of cosmic ray origin (e.g. Ptuskin 
(2001), see also discussion below). 

At the ejecta dominated stage which precedes the adi- 
abatic stage, the average spectrum of the run-away par- 
ticles is different from p~'^. Let us consider the general 
case and assume that p oc r~^, Rsh oc t^ and hence 



Pn 



sh 



■ , SO that 



i?m(p) oz p^-^-t: . Now Eq. (25) at A < 4/(6 - s) gives 
the following shape of the average spectrum of run-away 
particles: 



1 A(5-s)-2 

Q{p)ocp ^->-(^-i\ 



(27) 



Wsh OC t 



A-1 



The maximum momentum of accelerated 



A characteristic of the adiabatic regime is the relation 
A = and therefore Eq. (27) gives Q{p) oc p~'^ in 

agreement with Eq. (26). The Chevalier - Nadyozhin so- 
lution for the ejecta-dominated stage has A = j^. With 
the set of parameters excepted in Section 3, we have then 
Q{p) oc p~^-^ for the acceleration at the shock produced 
by the Type II SN in the presupernova star wind (s = 2, 
A; = 10, A = 7/8), and Q{p) oc p~'^ for the acceleration 
at the shock produced by the Type la SN in the uniform 
interstellar medium (s = 0, fc = 7, A = 4/7). Thus the cos- 
mic rays accelerated at the ejecta dominated stage have 
higher energies than at the later adiabatic stage but the 
average energy spectrum of produced cosmic rays is rather 
steep at the "canonical" choice of presupernova star pa- 
rameters. 

The results of our numerical calculations of the average 
spectra for Type II and Type la SNe are shown in Fig. 2. 
The parameter k is equal to 0.1 in a high- velocity regime 
(7) and it is equal to 0.04 in a low- velocity regime (9). 

The calculations for Type II SN are based on Eqs. (5), 
(6), (15), (17) and (25). For the set of parameters accepted 
in the present paper, the Type II supernovae are able to 
accelerate cosmic ray protons to the maximum energy of 
the order 4 x 10^^ eV if the acceleration starts one week 
after the SN explosion when Ush ~ 2.4 x 10^ km/s. The 
energy spectrum is close to p~^ at energies less than about 
6 X 10^^ eV and it experiences the steepening above this en- 
ergy. Thus the proton knee lies at about 6 x 10^^ eV in good 
agreement with the observational data. The sharp dip in 
the average proton spectrum at p/mc 1 x 10^ — 3 x 10^ 
is caused by the assumed abrupt change of the gas density 
at the boundary between the dense Red Super Giant wind 
and the low density bubble. We run calculations up to the 
maximum shock radius 60 pc (the corresponding SNR age 
is 9 x 10^ yr) when the Mach number approaches 3 and 
the use of the particle spectrum (1) characteristic of the 
strongly modified shocks can not be longer justified. The 
protons are accelerated to about 20 GeV at this moment. 

The calculations for Type la SN in Fig. 2 are based on 
Eqs. (5), (6), (14), (16) and (25). The calculations were 
made for the shock radius range from 0.2 to 30 pc (the 
SNR age from 4 yr to 1.3 x 10^ yr). The shock velocity 
is changing during this period from 2.8 x 10^ km/s to 
91 km/s. The protons arc accelerated from the maximum 
energy 7 x 10^^ eV to about 10 GeV. The approximate 
proton spectrum extends to the knee at about 3 x 10^^ 
eV. 

The average source spectrum produced by Type la SNe 
is multiplied by 1/4 in Fig. 2 that corresponds to their 
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Fig. 2. The solid line shows the average source spectrum Q{p)p'^c (given in units ^VgnS per steradian) for protons 
released into the interstellar medium during SNR evolution after SNII explosion in the wind of RSG progenitor star. 
The dashed line presents the case of SNIa explosion in the uniform interstellar gas; the average source spectrum is 
multiplied by 1/4. The dotted line shows the shape of proton source spectrum used by Horandel (2003) to fit the 
KASCADE data. 



relative burst rate and hence reflects the relative contri- 
bution of this type of supernovae to the total production 
of cosmic rays in the Galactic disk as compared to the 
contribution of Type II SNe. 

5. Consistency with Cosmic Ray Data and 
Discussion 

The spectrum of high-energy cosmic rays in the Galaxy is 
of the form / cx p~'^ , 7 = 7s + & under the steady state 
conditions when the action of cosmic ray sources (with 
the source power Q cx p ) is balanced by the escape of 
energetic particles from the Galaxy (with the escape time 
T cx p^^)- The observed at the Earth all-particle spectrum 
of cosmic rays is close to / cx p"**-^ at energies £' > 10 
Gev/nuclcon with a characteristic transition (the knee, 
Kulikov & Khristianscn 1958) ranging across less than one 
decade in the vicinity of 4 x 10^^ cV to the another power- 
law / cx p~^-^. The last extends to about 5x 10^^ eV where 
the second knee with the break ~ 0.3 is seen in the 



cosmic ray spectrum, see Horandel f2003') for review. This 
structure is usually associated with a severe decrease of the 
efficiency of cosmic ray acceleration or / and confinement in 
the Galaxy. The extragalactic component of cosmic rays 
probably dominates at > 3 x 10^^ eV (Gaisser et al. 
11993(1 . In the alternative interpretation (Berezinsky et al. 
12004(1 . the Galactic component falls steep (with 7 ~ 6) at 
E > 10^^ eV and the extragalactic component dominates 
from energy '--^ 3 x 10^^ eV and on. 

The exponent b ~ 0.3. ..0.7 was obtained from the 
data on the abundance of secondary nuclei at energies 
10^ to 10^^ eV/nucleon. The secondary nuclei are pro- 
duced in cosmic rays in a course of nuclear fragmentation 
of more heavy primary nuclei moving through the inter- 
stellar gas. The uncertainty in the value of b is mainly 
due to the choice of specific model of cosmic ray trans- 
port in the Galaxy, see Ptuskin 1(2001(1 . Hence it follows 
that the source exponent below the first knee lies in the 
range 7s = 4.0. ..4. 4. The value 7s « 4.0 for the aver- 
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age source spectrum was obtained above in the consider- 
ation of particle acceleration by SNR shocks during their 
adiabatic evolution (though smaller b ~ 0.3 and conse- 
quently larger 7^ ~ 4.4 would be more favorable for the 
explanation of high isotropy of cosmic rays observed at 
10^^ to 10^'' eV). According to the results of Section 4, 
the calculated average source spectrum for protons 
accelerated by a "typical" Type II SNe extends up to 
about 6 X lO-'^^eV that coincides with the observed po- 
sition of the knee '--^ 4 x 10^^ eV within the accuracy of 
our analysis. The knee position at 3-5 PeV was deter- 
mined in the recent KASCADE experiment (H. Ulrich et 
al. 12003(1 . The scaling of the knee position in our model is 
Pknee oc Z >c^£ M^^'^ M^^^ u^^ for the explosion in the stel- 
lar wind and Pknoc oc ZxS^SM^-^^^'^nQ^^ for the explosion 
in the uniform interstellar medium. 

As was reminded earlier, the diffusive shock accelera- 
tion at the strong nonmodified shock produces the spec- 
trum / cx p~^. The back reaction of efficiently accelerat- 
ing particles modifies the shock structure that results in 
a more flat particle spectrum (see references at the begin- 
ning of Section 2 and Eq. (1) where a ~ 0.5 if the shock 
modification is very strong). However, the numerical sim- 
ulation of acceleration by SNR shocks under the standard 
assumption of Bohm diffusion in the shock precursor (cal- 
culated for the interstellar magnetic field strength) and 
with efficient confinement of accelerated particles during 
all SNR evolution gives the overall source spectrum that 
is close to (Berezhko et al. I199()f) . Berezhko & Volk 
H2000|l pointed out that the last result is in some sense 
accidental. The late stages of the SNR evolution are im- 
portant here since relatively weak shock produces steep 
particle spectrum that has an effect on the overall spec- 
trum. The situation is different in the model discussed 
in the present paper because the coefficient of diffusion 
is strongly decreasing with SNR age and the cosmic rays 
with energies larger than 10 — 30 GeV/nucleon leave the 
supernova shell as the run away particles when the shock 
remains strong. The final average source spectrum of high- 
energy cosmic rays with energies larger than 10—30 GeV/n 
is close to provided that the shock evolution is approx- 
imately adiabatic and the efficiency of particle acceleration 
^cr is roughly constant. The source spectrum of particles 
with energies less than 10 — 30 Gev/n may be more steep 
because they are accelerated by not very strong shocks. In 
this connection it should be noted that the source spec- 
trum in the basic empirical model of cosmic ray propaga- 
tion in the Galaxy is of the form Q{p)p^ oc p^^ '* at S < 30 
GeV/n, and Q{p)p'^ cx p-2 i5 at _E > 30 GeV/n, see Jones 
et al. P)01,l . 

There are other important differences between the 
standard and the presented here models of cosmic ray 
acceleration. As noted before, our model of cosmic ray 
acceleration with strong increase in time of the diffusion 
coefficient and the corresponding decrease of maximum 
particle energy may naturally explain why the SNRs are 
generally not bright in very high energy gamma-rays at 



the age larger than a few thousand years. At this period 
of time, there arc no particles with energies needed to 
generate the very high energy gamma-rays in SNR shell. 
Another problem is the contribution of gamma-ray emis- 
sion from numerous unresolved SNRs with relatively flat 
spectra to the diffuse galactic background at very high en- 
ergies. Working in the standard model Berezhko & Volk 
( 200i|) took the maximum possible energy of protons ac- 
celerated in SNRs equal to 10^"* eV, i.e. well below the 
knee position, and it allowed not to exceed the upper lim- 
its on the diffuse gamma-ray emission at 4 x 10"'^^ to 10^"^ 
eV obtained in the Whipple, HEGRA, and TIBET exper- 
iments. In the model considered in the present work even 
with the efficient proton acceleration that may go beyond 
the knee, the expected gamma-ray emission from SNRs at 
4x 10^^ eV is order of magnitude smaller than in the model 
with Bohm diffusion. For a similar reason, the standard 
model compared to the present model predicts larger ratio 
of fluxes of secondary and primary nuclei formed at very 
high energies through the reacceleration of secondaries by 
strong shocks and through the direct production of sec- 
ondaries by primary nuclei with flat energy spectra inside 
SNRs, see Berezhko et al. C2Q03|. 

The interpretation of energy spectrum beyond the knee 
in the present model is associated with the cosmic ray 
acceleration during the ejecta dominated stage of SNR 
evolution when the protons gain by an order of magni- 
tude larger energy than at adiabatic stage but the num- 
ber of particles involved in the shock accelerated is rela- 
tively small. The average source spectrum of accelerated 
particles is not universal at this stage. It has a power law 
high-energy asymptotics with the exponent 7s which value 
is very sensitive to the parameter k. The last is not well 
determined from the observations but the typical values 
accepted in our calculations were A; = 10 for the Type II 
SN explosion in the wind of a Red Super Giant progenitor, 
and k = 7 for the Type la SN explosion in the uniform in- 
terstellar medium (see Chevalier & Liang IT^S^ that give 
7s = 6.5 and 7s = 7 respectively, see Section 4. To illus- 
trate the range of possible uncertainty, it is worth noting 
that the value k = 5.4 was suggested for the Type la SNe 
by Imshennik at al. (1981). This value of k results in 7s = 
4.3 at the ejecta dominated stage. 

The breaks and cutoffs in the spectra of ions with dif- 
ferent charges should occur at the same magnetic rigidity 
as for protons, i.e. at the same ratio p/Z (or E/Z for ul- 
trarelativistic nuclei) . The data of KASCADE experiment 
(Ulrich et al. "2003] for the most abundant groups of nuclei 
(protons, helium, CNO group, and the iron group nuclei) 
are, in general, consistent with this concept. According 
to Horandel (2003) the good fit to the observations is 
reached if an individual constituent ion spectrum has a 
gradual steepening by (57 ~ 2 at energy 4 x 10^^ Z eV. 
Eq (27) shows that the value S% = 2 can be obtained at 
A; = 9, s = 2 (the SN explosion in the progenitor wind), or 
A; = 6.6, 3 = (the SN explosion in the uniform interstel- 
lar medium) that is not very different from our accepted 
"typical" values, see Fig. 2. 
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At present, the main problem of the data interpre- 
tation centers around the second knee in the cosmic ray 
spectrum. The natural assumption that all individual ions 
has only one knee at ~ 4 x 10^^ Z eV and that the knee in 
the spectrum of iron {Z = 26) expected at about 10^^ eV 
explains the second knee in the all-particle spectrum does 
not agree with the observed position of the second knee 
at 5 X 10^^ eV. One way out was suggested by Horandel 
l|2()()3|l who included all elements up to Z = 92 into the 
consideration and assumed that 7s decreases with Z to 
rise the contribution of ultra heavy nuclei from Galactic 
sources to the cosmic ray flux at > lO^"^ eV. Of consider- 
able promise is the approach by Sveshnikova H2003fl who 
took into account the dispersion of parameters of SN ex- 
plosions in her calculations of the knee position and the 
maximum particle energy. It leads to the widening of the 
energy interval between the two knees in the overall all- 
particle spectrum. This analysis should be supplemented 
by the account of different chemical composition of the 
progenitor star winds that determines the composition 
of accelerated cosmic rays (Silberberg et al. I1991|l . We 
plan the future work on this topic in the frameworks of 
the model developed in the present paper and Paper I. It 
should be noted that the model by Bcrezinsky et al. 1)2004(1 
is quite consistent with the falling down of the flux from 
Galactic sources above 10^^ eV since the conjunction with 
the intergalactic cosmic ray flux in their model occurs at 
relatively low energy. 

There is also a very different scenario which assumes 
the strong reacceleration of cosmic rays above the knee by 
the collective effect of multiple SNR shocks in the violent 
regions of Galactic disk (Axford 119941 Bykov & Toptygin 
nUUT Klepach et al. IMlUjl or Galactic wind (Volk & 
Zirakashvih l^UMjl . 

Finally, it is worth noting that in principle the knee 
may arise not in the sources but in the process of cosmic 
ray propagation in the Galaxy, e.g. as a result of interplay 
between the ordinary and the Hall diffusion (Ptuskin et 
al. 119931 Roulet l^003|l . However, this explanation requires 
the existence of the power-law source spectrum which ex- 
tends without essential breaks up to about 10^^ eV or even 
further. 

6. Conclusion 

The accounting for non-linear effects which accompany 
the cosmic ray streaming instability raises the maximum 
energy of accelerated particles in young SNRs above the 
standard Bohm limit by about two orders of magnitude. 
It also considerably reduces the maximum energy of par- 
ticles that are present inside SNRs at the late Sedov stage 
if, as it was assumed in our calculations in Section 4, the 
cosmic ray diffusion coeflticient downstream of the shock is 
not much smaller than the diffusion cocfflcicnt in the cos- 
mic ray precursor of the shock and the energetic particles 
with p oc Pniax runs away from the SNR interior. In the 
present paper we studied the effect of arising strong time 
dependence of maximum particle momentum Pmax(i) on 



the average spectrum of cosmic rays injected into interstel- 
lar space from many supernova remnants over their life- 
time. The instantaneous cosmic ray spectrum at strongly 
modified shock is flat (/o cx p"'*"'"", a > 0, Eq. (1)) and 
the particle energy density is mainly determined by the 
particles with maximum momentum Pmax(i)- The instan- 
taneous source spectrum of the run away particles is close 
to the delta function {qi-a{t,p) oc (5(pi„ax(^) —p), Eq. (24)). 
At the same time, the assumption that the constant frac- 
tion of incoming gas momentum flux goes to the cosmic 
ray pressure at the shock, and the fact that the super- 
nova remnant evolution is adiabatic lead to the average 
on ensemble of SNRs source spectrum of ultrarelativistic 
particles that is close to Qra oc p^'* from energies 10 — 30 
GeV/n up to the knee position in the observed cosmic 
ray spectrum independent of the value of a, see Eq. (26) 
and Fig 2. This source spectrum is consistent with the 
empirical model of cosmic ray propagation in the Galaxy. 
The acceleration at the preceding ejecta-dominated stage 
of SNR evolution provides the steep power-law tail in the 
particle distribution at higher energies up to ~ 10^^ eV (if 
the iron nuclei dominate at these energies). The knee in 
the observed energy spectrum of cosmic rays at ~ 4x 10^^ 
eV is explained in our model by the transition from the 
particle acceleration at the ejecta-dominate stage to the 
adiabatic stage of SNR shock evolution. In spite of ap- 
proximate character of our consideration, it seems that 
the suggested scenario of particle acceleration can explain 
the energy spectrum of Galactic cosmic rays. 
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Appendix A: Tliin sliell approximation 

The thin-shell approximation can be used when the swept- 
up gas is concentrated in a thin layer behind the shock. 
In particular, it is applied to the case of a spherical adi- 
abatic shock, see Ostriker J. P. & McKee C.F. H1988;i and 
Bisnovatyi-Kogan & Silich H1995|l for detail. The total 
mass of the gas shell involved in the motion and confined 
by the shock of radius Rgh in the spherically symmetrical 



M = Mcj + 47r / drr^p{r), 



(A.l) 



where A/oj is the ejected mass, p is the density of ambient 
gas. 

The equation of momentum conservation is 



d{Mu) 
dt 



4^i?3\(Pi„-P) 



(A.2) 



Here u is the gas velocity behind the shock, Pjn is the 
pressure behind the shock, and P is the pressure of am- 
bient gas. For the adiabatic blastwave, u is related to the 
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shock velocity Ush = by the equation Ush = ^^H^' 



where 7ad is the ratio of the specific heats (adiabatic in- 
dex). The energy of explosion E = £th + \M%P' consists 

of the internal energy £th = 3(7^ ^--i) ^^"^ kinetic 
energy. 

Now for the case of very strong shock when can 
be omitted as negligible compared to P, Eq.(A2) can be 
presented as: 



d[Muf 12(7, 



ad 



1) 



rf-Rsh (7ad + 1)-Rsh 



EM 



\{Muf 



(A.3) 



The solution of Eq. (A3) allows finding the shock velocity 
and the shock age as functions of the shock radius: 



Ush(-Rsh) 

i(i?sh) = 

where w - 
main text. 



7ad 



1 



2 



-I 1/2 



drr'"-Hl{r) 







Ushir) ' 

6(7ad-l) 

7ad + l 



(A.4) 



that coincides with Eq. (17) in the 
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